clear all
path(path,'../TOOLS_01');

%%%%% le malha
%aux=load_gmsh4('malha_32_ele.msh',[1,2]);

aux=load_gmsh4('Malha_problema_direto_3_nos.msh',[1,2]);

n_nos=aux.nbNod;
n_elementos=aux.nbTriangles;
coords=aux.POS;
topologia_aux= aux.ELE_NODES(:,1:3);
n_eletrodos=max(aux.LINES(:,3));
topologia=topologia_aux(n_eletrodos+aux.nbLines+1:end,:);

ele_tags=aux.ELE_TAGS(n_eletrodos+aux.nbLines+1:end,:);
% idx de objetos
idx_objeto1=find(ele_tags(:,1)==2);
idx_objeto2=find(ele_tags(:,1)==3);
idx_objeto3=find(ele_tags(:,1)==4);

% idx de elementos
idx_malha=find(ele_tags(:,1)==1);

altura=0.05; % altura da malha
rho_eletrodo=0.5;
corriente=1e-3;

% rho_malha=zeros(1,n_elementos+n_eletrodos);
% rho_malha(1:n_elementos)=rho;
% rho_malha(n_elementos+1:end)=rho_eletrodo;

topologia_eletrodo=zeros(n_eletrodos,4);
for i=1:2:n_eletrodos*2
    topologia_eletrodo((i-1)/2+1,1:3)=[aux.LINES(i,1),aux.LINES(i,2),aux.LINES(i+1,2)];
    topologia_eletrodo((i-1)/2+1,4)=topologia_aux((i-1)/2+1);
end


%%%% objeto
rho=2.0;
rho_objeto=20.0;

rho_malha=zeros(1,n_elementos+n_eletrodos);
rho_malha(idx_malha)=rho;
rho_malha(idx_objeto2)=rho;
rho_malha(idx_objeto3)=rho;




rho_malha(idx_objeto1)=rho_objeto;
rho_malha(n_elementos+1:end)=rho_eletrodo;
V_objeto = p_direto(corriente,n_nos,n_elementos,coords,topologia,rho_malha,n_eletrodos,topologia_eletrodo,altura ,topologia_eletrodo);
%load lixo


%%%% valor inicial



rho=5.0;
rho_objeto=rho;


rho_malha=zeros(1,n_elementos+n_eletrodos);
rho_malha(idx_malha)=rho;
rho_malha(idx_objeto1)=rho;
rho_malha(idx_objeto2)=rho;
rho_malha(idx_objeto3)=rho;


rho_malha(n_elementos+1:end)=rho_eletrodo;
alpha=1e-13;
%%%% valor inicial


[ centroides ] = calc_centroides( coords, n_elementos,topologia);
max_x=max(centroides(:,1));
min_x=min(centroides(:,1));
max_y=max(centroides(:,2));
min_y=min(centroides(:,2));
res=64;
[XI,YI]=meshgrid(linspace(min_x,max_x,res),linspace(min_y,max_y,res));

%%%% matriz regurarizacao
norm_value_FtF=10.0;
[FtF]=calc_FtF_gaus(centroides,n_elementos,alpha,norm_value_FtF);



n_iteracoes=8;
passo_nr=1;

for i=1:n_iteracoes
    [Y_global,Y_local] = monta_rigidez_tri(n_nos,n_elementos,coords,topologia,rho_malha,n_eletrodos,topologia_eletrodo,altura);
    [ H,V_calculado ] = monta_jacobiano_2D(Y_global,Y_local,rho_malha,corriente,n_elementos,n_eletrodos,n_nos,coords,topologia_eletrodo);
    delta_rho=(H'*H+FtF)\(H'*(V_objeto-V_calculado)'-0.0*FtF*rho_malha(1:n_elementos)');
    rho_novo=rho_malha(1:n_elementos)+passo_nr*delta_rho';
    rho_novo(find(rho_novo<=0))=1.0;
    %%%%% plota imagem
     F = TriScatteredInterp(centroides(:,1),centroides(:,2),rho_novo');
     QZI=F(XI,YI);
     figure(2),surf(XI,YI,QZI)
     
     F1 = TriScatteredInterp(centroides(:,1),centroides(:,2),delta_rho);
     QZI1=F1(XI,YI);
     figure(3),surf(XI,YI,QZI1)
     
%      pause
    %%%%% plota imagem
    rho_malha(1:n_elementos)=rho_novo';
    
    erro(i)=norm((V_calculado-V_objeto))+rho_novo*FtF*rho_novo';
    if (i>4)
        if erro(i)>erro(i-3)
            passo_nr=passo_nr;
        end
    end
            
    figure(1),semilogy(erro,'.-'),title(sprintf('passo NR = %g',passo_nr))
end


F = TriScatteredInterp(centroides(:,1),centroides(:,2),rho_novo');
QZI=F(XI,YI);
figure(20),surf(XI,YI,QZI),colorbar


figure(1),semilogy(erro,'.-'),grid



